Análisis Accidentes de Tráfico NY

0.1 Cargar librerías y lectura de datos:

library(tidyverse)
library(viridis)
library(leaflet)
library(ggplot2)
library(ggbiplot)
library(dplyr)
library(GGally)
library(factoextra)
library(cluster)
library(tsibble)
library(psych)
library(caret)
library(MASS)
library(ca)
library(leaflet)
library(leaflet.extras)
library(stringdist)
library(igraph)
library(tidygraph)
library(ggraph)
library(pROC)
library(patchwork)
library(dendextend)
library(FactoMineR)
library(plotly)

#data <- read.csv("Vehicle_Collisions.csv")

0.2 Lectura datos (temporal)

Jon:

data <- read.csv("/DatosJon/UPNA/3º Ciencia de Datos/S6/Analisis Multivariante y Visualizacion de Datos/Trabajo/Datos/Vehicle_Collisions.csv")

Rubén:

#data <- read.csv("C:/Users/ruben/OneDrive/Escritorio/UPNA/CIENCIA DE DATOS/6ºSemestre/AnalisisMultivarianteVisualizacionDatos/Trabajo Final/Motor_Vehicle_Collisions_-_Crashes.csv")

Fermin:

# data <- read.csv("C:/Users/mintx/Mi unidad/CdD UPNA/.3.2/Analisis multivariante y visualizacion de datos/Trabajo/Vehicle_Collisions.csv")

1 INTRODUCCION

En nuestro conjunto de datos, se recoje información sobre accidentes de tráfico en la ciudad de Nueva York. La información proviene de la página de datos del gobierno de EEUU, específicamente del apartado correspondiente a la ciudad de Nueva York.

Cada individuo corresponde a un accidente de tráfico, que se representa como una fila en el conjunto de datos. Las columnas incluyen información como la fecha y hora del accidente, la ubicación, el tipo de vehículo, y los factores contribuyentes al accidente entre otros muchos. Los accidentes se han almacenado desde el 1/01/2013 hasta la actualidad.

Son accidentes en los que se ha rellenado un informe policial, por lo que podrían no incluirse algunos accidentes menores o accidentes que no han sido reportados a la policía. Ya que este informe es necesario para los accidentes donde hay fallecidos, o daños de al menos 1000$.

Además, como la forma de recoger los datos ha ido cambiando a lo largo de los años, es posible que hace años no se almacenase información como puede ser las coordenadas, y que en los últimos años con la implementación de dispositivos electrónicos se almacenen las coordenadas exactas de cada accidente. También habrá diferencias en la información dependiendo de la forma en la que el agente de policia correspondiente haya guardado los datos, ya que no son sensores sin fallo los que determinan la mayoría de variables (como pueden ser el tipo de vehículo o factor contribuyentes de la colisión).

Con todo esto, en nuestro trabajo trataremos de analizar las opciones que nos da el conjunto, e ir mostrando los resultados obtenidos visualmente. Y así poder encontrar patrones, agrupar los datos, o ver si podemos encontrar relación entre las variables que nos permitan entender mejor el problema tratado.

2 LIMPIEZA INICIAL Y MUESTREO

En primer lugar, observamos que el dataset contiene más de 2 millones de filas, pero que además muchas de las variables están prácticamente vacías o con información redundante. Además, dada la naturaleza de los datos geolocalizados, hay ejemplos que no contienen esa información y no nos son de ayuda.

Por tanto, vamos a realizar una primera limpieza del dataset, eliminando tanto variables redundantes o con poca información, como ejemplos que no nos aporten la información suficiente.

2.1 Eliminar variables:

• ZIP.CODE (ya tenemos la variable BOROUGH que contiene el barrio)
• ON.STREET.NAME, OFF.STREET.NAME, CROSS.STREET.NAME (es información que ya se almacena en las coordenadas)
• VEHICLE.TYPE.CODE.3, VEHICLE.TYPE.CODE.4, VEHICLE.TYPE.CODE.5 (prácticamente todos son valores nulos)
• CONTRIBUTING.FACTOR.VEHICLE.2, CONTRIBUTING.FACTOR.VEHICLE.3, CONTRIBUTING.FACTOR.VEHICLE.4, CONTRIBUTING.FACTOR.VEHICLE.5 (prácticamente todos son valores nulos)
• COLLISION_ID (es un identificador que no aporta información)
• LOCATION (es la combinación exacta de las variables LATITUDE Y LONGITUDE)

data_cleaned <- data |> 
  dplyr::select(-ZIP.CODE, -ON.STREET.NAME, -OFF.STREET.NAME, -CROSS.STREET.NAME, -VEHICLE.TYPE.CODE.3, -VEHICLE.TYPE.CODE.4, -VEHICLE.TYPE.CODE.5, -COLLISION_ID, -CONTRIBUTING.FACTOR.VEHICLE.5, -CONTRIBUTING.FACTOR.VEHICLE.4, -CONTRIBUTING.FACTOR.VEHICLE.3, -CONTRIBUTING.FACTOR.VEHICLE.2,-LOCATION)

Con el filtrado realizado, obtenemos un dataset mucho más manejable en cuanto a variables se refiere, que nos permitirá trabajar mejor el problema. Haciendo un pequeño resumen de las variables que nos quedan:

• CRASH.DATE: Fecha en la que se produció el accidente (Date)
• CRASH.TIME: Hora en la que se produció el accidente (Time)
• BOROUGH: Barrio de Nueva York donde sucedió el accidente (String)
• LATITUDE/LONGITUDE: Coordenadas geográficas (Int)
• NUMBER.OF.(PERSONS/PEDESTRIANS/CYCLIST/MOTORIST).(KILLED/INJURIED): Número y tipo de personas implicadas en el accidente, diferenciando heridos y fallecidos (Int)
• CONTRIBUTING.FACTOR.VEHICLE: Causa principal del accidente (String)
• VEHICLE.TYPE.CODE.1: Tipo del primero vehiculo implicado (String)
• VEHICLE.TYPE.CODE.2: Tipo del segundo vehiculo implicado (String)

2.2 Eliminar individuos:

• Individuos que no tengan LATITUDE y/o LONGITUDE (No nos sirven para localizar el accidente)
• Individuos que no tengan ningún Vehiculo implicado en VEHICLE.TYPE.CODE
• Individuos que no tengan barrio BOROUGH (ya que tenemos los suficientes que si tienen esa información)
• Individuos que no tengan una razon que haya causado el accidente CONTRIBUTING.FACTOR.VEHICLE
• Individuos con valor de localización erronea (Coordenadas 0) LONGITUDE

data_cleaned <- data_cleaned |> 
  filter(
    !is.na(LATITUDE) & !is.na(LONGITUDE),
    !(is.na(VEHICLE.TYPE.CODE.1) & is.na(VEHICLE.TYPE.CODE.2)),
    !(VEHICLE.TYPE.CODE.1 == "" & VEHICLE.TYPE.CODE.2 == ""),
    !(BOROUGH == ""),
    !(CONTRIBUTING.FACTOR.VEHICLE.1 == "Unspecified"),
    !(CONTRIBUTING.FACTOR.VEHICLE.1 == ""),
    !(LATITUDE==0),
    !(LONGITUDE==0))

2.3 Renombrar variables:

Antes de empezar a trabajar con los datos, vamos a renombrar ciertas variables para facilitar la programación.

data_cleaned <- data_cleaned |> 
  rename(
    DATE = CRASH.DATE,
    TIME = CRASH.TIME,
    NUM_PERSONS_INJURED = NUMBER.OF.PERSONS.INJURED,
    NUM_PERSONS_KILLED = NUMBER.OF.PERSONS.KILLED,
    NUM_PEDESTRIANS_INJURED = NUMBER.OF.PEDESTRIANS.INJURED,
    NUM_PEDESTRIANS_KILLED = NUMBER.OF.PEDESTRIANS.KILLED,
    NUM_CYCLIST_INJURED = NUMBER.OF.CYCLIST.INJURED,
    NUM_CYCLIST_KILLED = NUMBER.OF.CYCLIST.KILLED,
    NUM_MOTORIST_INJURED = NUMBER.OF.MOTORIST.INJURED,
    NUM_MOTORIST_KILLED = NUMBER.OF.MOTORIST.KILLED,
    CAUSE = CONTRIBUTING.FACTOR.VEHICLE.1,
    VEHICLE_1 = VEHICLE.TYPE.CODE.1,
    VEHICLE_2 = VEHICLE.TYPE.CODE.2
  )

2.4 Eliminación de filas incompletas:

2.4.1 Personas afectadas:

Vamos a comprobar si el número de personas heridas o fallecidas coincide con el número de heridos o fallecidos que se han dado en la fila. Para ello, vamos a crear una tabla con los diferentes tipos de heridos y fallecidos, y luego vamos a comprobar si el número total coincide con el número de personas heridas o fallecidas que se han dado en la fila.

type_injured <- data_cleaned |> 
  dplyr::select(NUM_PERSONS_INJURED, NUM_PEDESTRIANS_INJURED, NUM_CYCLIST_INJURED, NUM_MOTORIST_INJURED,
         NUM_PERSONS_KILLED, NUM_PEDESTRIANS_KILLED, NUM_CYCLIST_KILLED, NUM_MOTORIST_KILLED) |> 
  pivot_longer(cols = everything(), names_to = 'TYPE_OF_INJURED', values_to = 'NUMBER')

# Calculamos el total de cada tipo:
type_injured |> 
  group_by(TYPE_OF_INJURED) |> 
  summarise(NUMBER = sum(NUMBER)) |> 
  arrange(desc(NUMBER))

# Ahora mostramos los indices de los que no coinciden:
filas_erroneas <- data_cleaned |> 
  dplyr::select(NUM_PERSONS_INJURED, NUM_PEDESTRIANS_INJURED, NUM_CYCLIST_INJURED, NUM_MOTORIST_INJURED,
         NUM_PERSONS_KILLED, NUM_PEDESTRIANS_KILLED, NUM_CYCLIST_KILLED, NUM_MOTORIST_KILLED) |> 
  mutate(ROW_ID = row_number(),
         TOTAL_NUMBER_INJURED = rowSums(across(c(NUM_PEDESTRIANS_INJURED, NUM_CYCLIST_INJURED, NUM_MOTORIST_INJURED))),
         TOTAL_NUMBER_KILLED = rowSums(across(c(NUM_PEDESTRIANS_KILLED, NUM_CYCLIST_KILLED, NUM_MOTORIST_KILLED)))) |>
  dplyr::select(ROW_ID, NUM_PERSONS_INJURED, TOTAL_NUMBER_INJURED, NUM_PERSONS_KILLED, TOTAL_NUMBER_KILLED) |> 
  filter(NUM_PERSONS_INJURED != TOTAL_NUMBER_INJURED | NUM_PERSONS_KILLED != TOTAL_NUMBER_KILLED)

head(filas_erroneas)

Podemos ver como existen filas que no cumplen, por ello eliminamos las filas erroneas, aquellas filas cuyo numero de personas heridas o fallecidas no coincide con el total de heridos o fallecidos que se han dado en la fila.

indices_erroneos <- filas_erroneas |> 
  pull(ROW_ID)

data_cleaned <- data_cleaned[-indices_erroneos, ]

2.5 Eliminacion de valores incorrectos de la variable ‘CAUSE’:

data_cleaned <- data_cleaned |> 
  mutate(CAUSE = recode(CAUSE,
    "Cell Phone (hand-Held)" = "Cell Phone (hand-held)",
    "Drugs (Illegal)" = "Drugs (illegal)",
    "Illnes" = "Illness",
    "Other Electronic Device" = "Other Devices",
    "Drugs (illegal)" = "Drugs (illegal)",
    "Unspecified" = "Unspecified",
    "Pedestrian/Bicyclist/Other Pedestrian Error/Confusion" = "Pedestrian Error/Confusion",
    "Aggressive Driving/Road Rage" = "Aggressive Driving"
  )) |> 
  filter(CAUSE != 1, CAUSE != 80) 

2.6 Formatear variable ‘DATE’:

data_cleaned <- data_cleaned |> 
  mutate(DATE = as.Date(DATE, format = "%m/%d/%Y"))

2.7 Estandarizar strings de las variables:

Para facilitar el trabajo con los datos y evitar valores duplicados, vamos a estandarizar las variables CAUSE, VEHICLE_1 y VEHICLE_2 convirtiéndolas completamente a mayúsculas.

data_cleaned$VEHICLE_1 <- toupper(data_cleaned$VEHICLE_1)
data_cleaned$VEHICLE_2 <- toupper(data_cleaned$VEHICLE_2)
data_cleaned$CAUSE <- toupper(data_cleaned$CAUSE)

2.8 Muestreo:

Una vez eliminados los individuos que no son de interes, el dataset se reduce a 888068 filas. Seguimos teniendo demasiadas observaciones, por lo que para trabajar con los datos de forma más comoda hemos decidido hacer un muestreo de 20.000 observaciones. Cabe destacar que si bien para tener los mismos resultados fijamos una semilla, durante el trabajo en ciertos analisis hemos tomado muestras diferentes para asegurarnos de que los resultados eran similares y no había demasiada variación debido a la muestra.

set.seed(123)
data_sampled <- sample_n(data_cleaned, 20000)

Como a partir de ahora trabajaremos con esta muestra en lugar del dataset completo, por comodidad lo almacenaremos en un nuevo archivo para solo tener que leer este.

write.csv(data_sampled, "data_sampled.csv", row.names = FALSE)
write.csv(data_sampled, "App/data_sampled.csv", row.names = FALSE)

2.8.1 Leer muestra:

data_sampled <- read.csv("data_sampled.csv")

3 ANALISIS DESCRIPTIVO

Una vez obtenido el dataset con el que trabajaremos, vamos a visualizar las diferentes variables y sus comportamientos para tratar de entender mejor con que datos estamos trabajando antes de aplicar las diferentes técnicas.

Además, también mostraremos la matriz de correlación entre las variables numéricas, para ver si hay alguna relación entre ellas.

Antes de empezar con las visualizaciones, vamos a crear un tema personalizado para los gráficos que vamos a mostrar. De esta manera, todos los gráficos tendrán el mismo estilo y será más fácil de leer.

3.1 Tema personalizado

Tanto en el informe como en la app, hemos definido un tema personalizado con colores característicos de Nueva York y particularmente los vehiculos, en el que destaca la presencia del amarillo color Taxi, y los grises del asfalto para dar contraste.

theme_nyc <- function() {
  theme_minimal(base_family = "Roboto Condensed") +
    theme(
      plot.background = element_rect(fill = "#F4F4F4", color = NA),
      panel.background = element_rect(fill = "#F4F4F4", color = NA),
      panel.grid.major = element_line(color = "#DADADA"),
      panel.grid.minor = element_blank(),
      axis.title = element_text(color = "#2E2E2E", size = 16, face = "bold"),
      axis.text = element_text(color = "#2E2E2E"),
      plot.title = element_text(color = "#1F3B73", size = 20, face = "bold", hjust = 0.5),
      plot.subtitle = element_text(color = "#2E2E2E", hjust = 0.5),
      legend.background = element_rect(fill = "#F4F4F4"),
      legend.key = element_rect(fill = "#F4F4F4"),
      legend.text = element_text(color = "#2E2E2E"),
      strip.background = element_rect(fill = "#FFD100"),
      strip.text = element_text(color = "#2E2E2E", face = "bold")
    )
}

3.2 Variable ‘DATE’:

Vamos a hacer una representación dual, en la que combinamos el gráfico de barras y el gráfico de líneas, para mostrar la frecuencia de accidentes según el mes o año. Para ello, vamos a crear una nueva variable que contenga la fecha en formato mes y otra en año.

data_time <- data_sampled |>
  mutate(DATE_MONTH = yearmonth(DATE),
         DATE_YEAR = year(DATE)) |>
  add_count(DATE, name = "FREQ_DAY") |>
  add_count(DATE_MONTH, name = "FREQ_MONTH") |>
  add_count(DATE_YEAR, name = "FREQ_YEAR") |>
  dplyr::select(DATE, DATE_MONTH, DATE_YEAR, FREQ_DAY, FREQ_MONTH, FREQ_YEAR, BOROUGH) |>
  arrange(DATE)


# Frecuencia de datos mensuales:
data_time |> 
  group_by(DATE_MONTH) |>
  summarise(FREQ = sum(FREQ_MONTH, na.rm = TRUE), .groups = "drop") |>
    ggplot(aes(x = DATE_MONTH, y = FREQ)) +
    geom_bar(stat = "identity", fill = "#1F3B73", alpha = 0.8) +
    geom_line(color = "#1F3B73", size = 1.2) +
    geom_point(color = "#FF4C4C", size = 2) +
    labs(title = "Frecuencia de Accidentes por Mes", x = "Fecha", y = "Frecuencia") +
    theme_nyc()

# Frecuencia de datos anuales:
data_time |> 
  group_by(DATE_YEAR) |>
  summarise(FREQ = sum(FREQ_YEAR, na.rm = TRUE), .groups = "drop") |>
    ggplot(aes(x = DATE_YEAR, y = FREQ)) +
    geom_bar(stat = "identity", fill = "#1F3B73", alpha = 0.8) +
    geom_line(color = "#1F3B73", size = 1.2) +
    geom_point(color = "#FF4C4C", size = 2) +
    labs(title = "Frecuencia de Accidentes por Año", x = "Fecha", y = "Frecuencia") +
    theme_nyc()

En estos gráficos observamos como ha ido evolucionando la frecuencia de accidentes a lo largo del tiempo, ya sea con datos acumulados mensuales o anuales. Se puede ver claramente como en el año del covid la frecuencia de accidentes disminuye significativamente, debido a que la gente no se movía tanto como en años anteriores. Sin embargo, la frecuencia de accidentes se mantiene casi constante, incluso llegando a disminuir en en 2022. Esto tambien se puede deber a que a raiz del covid la gente ha cambiado su forma de moverse, y ahora se mueven más en bicicleta o andando, lo que puede hacer que haya más accidentes de este tipo. También, hemos comprobado que después del covid, en EEUU, modificaron leyes de tráfico que hizo que la gente sea mas consciente de la seguridad vial, y por lo tanto, haya menos accidentes.

A continuación, vamos a representar la frecuencia total acumulada por cada mes. Es decir, la suma de frecuencia para cada uno de los meses en todos los años pasados:

data_time |> 
  group_by(DATE_MONTH) |> 
  summarise(FREQ_MONTH = sum(FREQ_MONTH)) |> 
  mutate(MONTH = month(DATE_MONTH, label = TRUE)) |> 
  ggplot(aes(x = MONTH, y = FREQ_MONTH)) + 
    geom_histogram(stat = "identity", fill = "#1F3B73", alpha = 0.8) +
    labs(title = "Frecuencia de accidentes acumulada por mes", x = "Fecha", y = "Frecuencia") +
    theme_nyc()

En este gráfico observamos cuales son los meses que máyor frecuencia de accidentes tienen. Observamos como a principios de verano es cuando más alta es la frecuencia (mayo, junio y julio) tal vez porque la gente se mueve más por el buen tiempo, pero aún teniendo que ir a trabajar. En agosto y septiembre se reduce la frecuencia, tal vez por el efecto de las vacaciones. Por otra parte, sin duda los meses de febrero (por tener 29 días) y abril (tal vez por ser el mes en el que más afecto el confinamiento de 2020) son los que menor frecuencia de accidentes tienen.

3.3 Variable ‘BOROUGH’:

En el siguiente apartado, vamos a tratar de ver la distribución de los accidentes por distrito, y su evolución con el paso del tiempo:

data_sampled |> 
  ggplot(aes(x = BOROUGH)) +
  geom_bar(fill = "#1F3B73") +
  geom_text(aes(label = after_stat(count)), stat = "count", vjust = -0.5, size = 3) +
  labs(title = "Distribución de Accidentes por Distrito",
       x = "Barrio", y = "Número de Accidentes") +
  theme_nyc()

Observamos como Brooklyn Queens y Manhattan acumulan el mayor número de accidentes. Esto concuerda con que son los distritos más grandes y céntricos de Nueva York. Se observa como Bronx y State Island muestran una frecuencia mucho más inferior, ya que son distritos periféricos.

Vamos a ver la evolución de la frecuencia de accidentes por cada distritos:

data_time |>
  group_by(DATE_MONTH, BOROUGH) |>
  summarise(accidents = sum(FREQ_DAY, na.rm = TRUE), .groups = "drop") |> 
  ggplot(aes(x = DATE_MONTH, y = accidents, color = BOROUGH)) +
    geom_line(size = 1.1) +
    geom_point(size = 1.5) +
    labs(title = "Evolución Temporal de Accidentes por Distrito",
         x = "Fecha", y = "Número de Accidentes", color = "Distrito") +
    theme_nyc() +
    scale_color_brewer(palette = "Set1") +
    theme(legend.position = "bottom")

Aqui observamos el mismo patrón de frecuencia que hemos analizado anteriormente, donde baja drásticamente a partir del covid. Además, parece que el comportamiento de todos los distritos es similar, y aumentan o disminuyen de manera correlacionada. También parece que el distrito del Bronx no disminuye tanto la frecuencia a partir de 2020, mostrándose por encima de Manhattan varias veces. También podemos observar como Manhattan era el distrito con mayor frecuencia hasta 2016, y desde 2017 hasta 2020 cede ese puesto a Brooklyn y Queens.

Por otrol lado, también podemos observar como Staten Island es el barrio que menos accidentes tiene, y además es el que menos ha disminuido la frecuencia a partir de 2020. Esto puede ser debido a que es un barrio más residencial, con menos tráfico y con menos pobalción, por lo que los accidentes son menos frecuentes.

3.4 Variables ‘LATITUDE’ y ‘LONGITUDE’:

En cuanto a las variables de localización, al ser un dataset con datos geográficos uno de los análisis importantes que debemos hacer es situando los accidentes en su ubicación correspondiente. Esto nos permite ver como están distrubuidos en la ciudad los accidentes.

En primer lugar, representaremos cada accidente mediante un punto rojo:

Observamos como prácticamente la totalidad del mapa queda cubierta de punto rojos, si bien es cierto que en distritos como Staten Island se aprecia menos densidad, tal y como intuíamos de los histogramas realizados previamente.

Para tratar de ver algo más claro, vamos a representar un mapa de calor, en el que las zonas más calientes representan aquellos accidentes en los que ha habido heridos o fallecidos (que ponderan de forma más significativa):

leaflet(data_sampled) %>%
  addTiles() %>%
  setView(lng = -73.9, lat = 40.7128, zoom = 10) %>%
  addHeatmap(
    lng = ~LONGITUDE,
    lat = ~LATITUDE,
    intensity = ~ 5*NUM_PERSONS_KILLED + NUM_PERSONS_INJURED,
    radius = 15,
    blur = 20,
    max = 0.05,
    group = "heat"
  )

Este mapa nos ofrece información más relevante. Aunque Manhattan concentra un gran número de accidentes, no presenta una alta densidad de heridos o fallecidos, como cabría esperar. Esto podría deberse a que la mayoría de los siniestros ocurren a baja velocidad en zonas urbanas. En contraste, distritos como Staten Island muestran una mayor densidad de accidentes con heridos, a pesar de registrar una menor frecuencia total de accidentes.

En la aplicación hemos desarrollado los mapas de forma interactiva, pudiendo aplicar diferentes filtros e intervalos de tiempo variables, para poder ver de una forma más clara y dinámica la evolución de los accidentes a lo largo del tiempo.

3.5 Variables de Conteo de Personas Afectadas:

Ahora, podemos comprobar el número de peatones, ciclistas y motoristas que han resultado heridos o muertos en los accidentes, junto al total de heridos y muertos:

    datos_resumen <- data_sampled |> 
      summarise(
        Total_Heridos = sum(NUM_PERSONS_INJURED, na.rm = TRUE),
        Heridos_Pedestrians = sum(NUM_PEDESTRIANS_INJURED, na.rm = TRUE),
        Heridos_Cyclists = sum(NUM_CYCLIST_INJURED, na.rm = TRUE),
        Heridos_Motorists = sum(NUM_MOTORIST_INJURED, na.rm = TRUE),
        Total_Muertos = sum(NUM_PERSONS_KILLED, na.rm = TRUE),
        Muertos_Pedestrians = sum(NUM_PEDESTRIANS_KILLED, na.rm = TRUE),
        Muertos_Cyclists = sum(NUM_CYCLIST_KILLED, na.rm = TRUE),
        Muertos_Motorists = sum(NUM_MOTORIST_KILLED, na.rm = TRUE)
      )
    
    datos_long <- datos_resumen |> 
      pivot_longer(cols = everything(),
                   names_to = c("Estado", "Tipo"),
                   names_sep = "_",
                   values_to = "Total")
    
    ggplot(datos_long, aes(x = Tipo, y = Total, fill = Tipo)) +
      geom_col(show.legend = TRUE) +
      geom_text(aes(label = Total), vjust = -0.5, size = 4, color = "#2E2E2E") +
      facet_wrap(~Estado, scales = "free_y") +
      labs(
        title = "Total de personas heridas y muertas por tipo",
        x = NULL,
        y = "Total",
        fill = "Tipo de persona"
      ) +
      theme_nyc() +
      theme(
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "bottom"
      )

En este gráfico observamos que en general, la frecuencia de heridos es mucho mayor respecto a fallecidos (6472 frente a 23). También observamos que el grupo más vulnerable son los motoristas, tanto en cuanto a heridos como a fallecidos se refiere. Seguidos de los peatones y ciclistas. Si bien, considerando que el número de ciclistas es muy inferior al de peatones, y apenas son la mitad de heridos, tal vez podríamos decir que ser ciclista si que conlleva un mayor riesgo de resultar herido en proporción.

Además, pese a ser variables enteras, como son numéricas, puede ser interesante mostrar un histograma de cada una de ellas para ver si hay normalidad en los datos:

vars <- c(
  "NUM_PERSONS_INJURED", "NUM_PEDESTRIANS_INJURED", "NUM_CYCLIST_INJURED", "NUM_MOTORIST_INJURED",
  "NUM_PERSONS_KILLED", "NUM_PEDESTRIANS_KILLED", "NUM_CYCLIST_KILLED", "NUM_MOTORIST_KILLED"
)

hist_list <- lapply(vars, function(var) {
  ggplot(data_sampled, aes_string(x = var)) +
    geom_histogram(bins = 30, fill = "#1F77B4", color = "white") +
    labs(title = var, x = NULL, y = NULL) +
    theme_nyc()
})

wrap_plots(hist_list, ncol = 4)

Observamos como la distribución en todas las variables llega al máximo en el 0, y va reduciéndose drásticamente. Lo cual es comprensible ya que se trata del número de personas heridas en accidentes, y en los accidentes leves no suele haber ningún herido.

Esto nos indica que no tenemos normalidad, lo cual habrá que tenerlo en cuenta a la hora de aplicar diferentes técnicas de análisis multivariante.

En el siguiente apartado, veremos la matriz de correlaciones entre todas las variables de conteo de personas afectadas. Como no tenemos normalidad, aplicaremos la matriz de correlaciones de Spearman:

data_short <- data_sampled[, c("NUM_PERSONS_INJURED", "NUM_PEDESTRIANS_INJURED", "NUM_CYCLIST_INJURED", "NUM_MOTORIST_INJURED", "NUM_PERSONS_KILLED",   "NUM_PEDESTRIANS_KILLED", "NUM_CYCLIST_KILLED", "NUM_MOTORIST_KILLED")]

colnames(data_short) <- c("Total Her.", "Peatones Her.", "Ciclistas Her.", "Motoristas Her.", "Total Fall.", "Peatones Fall.", "Ciclistas Fall.", "Motoristas Fall.")

ggpairs(
  data_short,
  title = "Matriz de correlaciones Spearman",
  upper = list(continuous = wrap("cor", size = 4, method = "spearman")),
  lower = list(continuous = wrap("points", alpha = 0.5)),
  diag = list(continuous = wrap("densityDiag"))
) +
  theme_nyc()

En conclusión, la matriz de correlaciones solamente muestra relaciones fuertes entre variables agregadas y sus subcomponentes, como entre personas heridas y motoristas heridos (ya que dentro de personas heridas están los motoristas, y además hemos visto que es el tipo que más contribuye). Sin embargo esta correlación es evidente, y no aporta ninguna información de valor.

En cuanto a las variables por separado sin datos repetidos, ninguna muestra correlación significativa con las demás, lo cual sugiere que no se comportan del mismo modo los tipos diferentes de implicados, y que no podemos explicar el número de heridos/fallecidos de ningún tipo basandónos en otro de los tipos.

3.6 Variable ‘CAUSE’:

Primero, vamos a ver la frecuencia de las distintas causas de los accidentes registrados. Para ello, vamos a crear una tabla con los diferentes factores contribuyentes y su frecuencia:

causes_cleaned <- data_sampled |> 
  count(CAUSE, name = "FREQUENCY") |>  
  arrange(desc(FREQUENCY))

Nos ayudamos de un gráfico de barras para mostrar las frecuencias.

causes_cleaned |> 
  ggplot(aes(x = reorder(CAUSE, FREQUENCY), y = FREQUENCY)) + 
  geom_bar(stat = "identity", fill = "steelblue") + 
  coord_flip() +
  labs(title = "Frecuencia de Factores Contribuyentes de los Accidentes",
       x = "Causas del Accidente",
       y = "Frecuencia") +
  theme_nyc()

Se observa que una gran parte de los accidentes se debe a distracciones. Además, aunque existen otras causas con una frecuencia significativa, muchas de las causas menos comunes apenas concentran accidentes. Esto sugiere la posibilidad de aplicar un análisis de agrupamiento (clustering) en el futuro, con el objetivo de reducir y simplificar el número de causas en grupos más representativos.

3.7 Variables ‘VEHICLE_1’ y ‘VEHICLE_2’:

Por ultimo, analizaremos las variables correspondientes al tipo de vehículo principal y secundario implicados en el accidente. De forma similar a la variable anterior, es una variable categórica, por lo que primero vamos a crear una tabla para cada variable con su frecuencia:

V1_cleaned <- data_sampled |> 
  count(VEHICLE_1, name = "FREQUENCY_V1") |>  
  arrange(desc(FREQUENCY_V1))

V2_cleaned <- data_sampled |> 
  count(VEHICLE_2, name = "FREQUENCY_V2") |>  
  arrange(desc(FREQUENCY_V2))

Ahora vamos a graficar los tipos de vehículos y sus frecuencias para VEHICLE_1 y VEHICLE_2. Para ello, vamos a usar un gráfico de barras. Comenzamos con la primera variable:

V1_cleaned |> 
  filter(FREQUENCY_V1 > 5) |> 
  ggplot(aes(x = reorder(VEHICLE_1, FREQUENCY_V1), y = FREQUENCY_V1)) + 
  geom_bar(stat = "identity", fill = "steelblue") + 
  coord_flip() +
  labs(title = "Frecuencia de Tipos de Vehículos (Primer involucrado)",
       x = "Vehículos",
       y = "Frecuencia") +
  theme_nyc()

Hemos mostrado solamente los vehiculos que tienen una frecuencia de almenos 5, ya que los inferiores generalmente se tratan de errores de tipografía, y por no extender en exceso el gráfico.

De manera similar, mostramos la segunda variable:

V2_cleaned |> 
  filter(FREQUENCY_V2 > 5, VEHICLE_2 != '') |> 
  ggplot(aes(x = reorder(VEHICLE_2, FREQUENCY_V2), y = FREQUENCY_V2)) + 
  geom_bar(stat = "identity", fill = "steelblue") + 
  coord_flip() +  # Rota el gráfico para mejor visualización
  labs(title = "Frecuencia de Tipos de Vehículos (Primer involucrado)",
       x = "Vehículos",
       y = "Frecuencia") +
  theme_nyc()

En ambos casos observamos que el vehículo más implicado son los coches Sedan habituales como cabía esperar, y va disminuyendo el número drásticamente al resto de tipos de vehículos.

4 METODOLOGÍA MULTIVARIANTE

4.1 PCA

Vamos a realizar un análisis de componentes principales (PCA) para reducir la dimensionalidad de los datos y ver si podemos encontrar patrones en ellos. Sin embargo, lo primero que vamos hacer es comprobar si tiene sentido hacer un PCA con los datos numéricos que tenemos, ya que no existen demasiadas variables y ya hemos podido comprobar que no presentan normalidad. Además, en la matriz de correlación ya hemos visto que no hay correlaciones significativas en los datos, lo cual puede complicarnos realizar las combinaciones necesarias para reducir el número de componentes principales.

Sin embargo, para tomar la decisión con un criterio más preciso, hemos decidido aplicar las metodologías de adecuación de la muestra que vimos para el análisis factorial mediante el test de KMO. Nos puede dar una idea también al trabajar con PCA, ya que también necesitamos correlación entre las variables, y explicar una parte significativa de la varianza en ambos métodos.

Aplicamos el test de KMO sobre las variables numéricas. Un valor menor o igual a 0.5 indica que los datos no son adecuados:

datos_numeric <- data_sampled |> 
  dplyr::select(NUM_PEDESTRIANS_INJURED, NUM_CYCLIST_INJURED, NUM_MOTORIST_INJURED,
         NUM_PEDESTRIANS_KILLED, NUM_CYCLIST_KILLED, NUM_MOTORIST_KILLED)

KMO(datos_numeric)

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = datos_numeric)
## Overall MSA =  0.47
## MSA for each item = 
## NUM_PEDESTRIANS_INJURED     NUM_CYCLIST_INJURED    NUM_MOTORIST_INJURED 
##                    0.47                    0.47                    0.48 
##  NUM_PEDESTRIANS_KILLED      NUM_CYCLIST_KILLED     NUM_MOTORIST_KILLED 
##                    0.44                    0.44                    0.53

Observamos el valor MSA general es de 0.47, por lo que no parece interesante realizar un PCA tal y como intuíamos. Sin embargo, vamos a realizar un gráfico de PCA sobre las dos componentes principales, para ver si aun así podemos extraer algún comportamiento interesante.

datos_pca <- data_sampled |> 
  select(NUM_PEDESTRIANS_INJURED, NUM_CYCLIST_INJURED, NUM_MOTORIST_INJURED,
         NUM_PEDESTRIANS_KILLED, NUM_CYCLIST_KILLED, NUM_MOTORIST_KILLED)

pca_result <- prcomp(datos_pca, center = TRUE, scale. = TRUE)

fviz_pca_var(pca_result, repel = TRUE, labelsize = 3, col.var = "contrib") +
  scale_color_viridis(option = "C") +
  theme_nyc()

Observamos que usando dos componentes, el porcentaje de variabilidad explicada apenas supera el 35%, lo cual no es lo ideal. Podemos extraer que las variables de heridos contribuyen de una manera más significativa que las variables sobre fallecidos. También podemos observar que el número de motoristar heridos y fallecidos presentan la misma dirección, lo que propone que, a mayor numero de heridos, mayor número de fallecidos. En cuanto a las otras variables, observamos que toman la dirección más opuesta posible, lo que nos indica que la correlación entre ellas es la mínima posible. Estas conclusiones respaldan el hecho de que la variabilidad explicada por las componentes principales sea tan baja.

Vamos a mostrar el gráfico de la varianza explicada por cada componente, para ver como va creciendo hasta el 100%, que se dará con las 6 variables:

pca_summary <- summary(pca_result)

var_exp_cum <- pca_summary$importance[3, ] * 100

data.frame(Componente = 1:length(var_exp_cum),
           Varianza_Acumulada = var_exp_cum) %>%
  ggplot(aes(x = Componente, y = Varianza_Acumulada)) +
  geom_bar(stat = "identity", fill = "steelblue", alpha = 0.6) +
  geom_line(color = "red", size = 1) +
  geom_point(color = "red", size = 3) +
  labs(title = "Varianza acumulada explicada por cada componente",
       x = "Componente Principal",
       y = "Varianza acumulada (%)") +
  theme_nyc()

Podemos observar que prácticamente todas las componentes aportan la misma cantidad de variabilidad explicada, debido a que no hay correlacion entre las variables como hemos observado antes. Además, observamos que con 5 componentes apenas llegamos al 85% de variabilidad explicada, por lo que una vez más no parece que podamos reducir la dimensionalidad de forma efectiva mediante PCA, ya que perderíamos demasiada información.

Por último, vamos a comprobar si es interesante analizar la posibilidad de realizar PCA separando las variables por heridos y muertos en lugar de todos juntos. Realizaremos el test de KMO del mismo modo que ahora, para ver si en estos casos sería interesante con dicha muestra.

4.1.1 PCA con heridos

datos_heridos <- data_sampled |> 
  dplyr::select(NUM_PEDESTRIANS_INJURED, NUM_CYCLIST_INJURED, NUM_MOTORIST_INJURED)

datos_heridos_scaled <- scale(datos_heridos)

KMO(datos_heridos_scaled)

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = datos_heridos_scaled)
## Overall MSA =  0.47
## MSA for each item = 
## NUM_PEDESTRIANS_INJURED     NUM_CYCLIST_INJURED    NUM_MOTORIST_INJURED 
##                    0.47                    0.47                    0.48

De manera similar al caso general, obtenemos que no es interesante realizar un PCA, ya que el KMO vuelve a rechazar la idea de realizar el PCA, con el mismo valor de 0.47, inferior a 0.5.

4.1.2 PCA con muertos

datos_heridos <- data_sampled |> 
  dplyr::select(NUM_PEDESTRIANS_KILLED, NUM_CYCLIST_KILLED, NUM_MOTORIST_KILLED)

datos_heridos_scaled <- scale(datos_heridos)

KMO(datos_heridos_scaled)

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = datos_heridos_scaled)
## Overall MSA =  0.5
## MSA for each item = 
## NUM_PEDESTRIANS_KILLED     NUM_CYCLIST_KILLED    NUM_MOTORIST_KILLED 
##                    0.5                    0.5                    0.5

Y por último, tampoco es interesante con los datos de solamente los fallecidos, obteniendose un KMO de 0.5. Como comentábamos al principio, son resultados esperables al haber observado previamente la matriz de correlaciones y la naturaleza propia de nuestras variables numéricas, por lo que deberemos centrarnos en realizar otro tipo de técnicas más apropiadas a nuestro problema.

Es interesante mencionar que, teniendo en cuenta los valores obtenidos en los tests de KMO, no vamos a considerar tampoco la posibilidad de realizar un análisis factorial, ya que el resultado del test es un claro indicador de que no es apropiado.

4.2 Clusterizaciones

4.2.1 Análisis Cluster Jerárquico en variable ‘BOROUGH’:

4.2.1.1 1. Seleccionar variables numéricas relevantes:

data_cluster <- data_sampled |> 
  dplyr::select(NUM_PERSONS_INJURED, NUM_PEDESTRIANS_INJURED, NUM_CYCLIST_INJURED, NUM_MOTORIST_INJURED,
         NUM_PERSONS_KILLED, NUM_PEDESTRIANS_KILLED, NUM_CYCLIST_KILLED, NUM_MOTORIST_KILLED)

# Escalar (normalizar) los datos
data_scaled <- scale(data_cluster)

4.2.1.2 2. Determinar el número óptimo de clusters:

wss <- vector()
set.seed(123)
for (k in 1:10) {
  km <- kmeans(data_scaled, centers = k, nstart = 5)  # menos repeticiones
  wss[k] <- km$tot.withinss
}

plot(1:10, wss, type = "b", pch = 19,
     xlab = "Número de clusters",
     ylab = "Suma de cuadrados intra-clúster (WSS)",
     main = "Método del codo (optimizado)")

4.2.1.3 3. Realizar el clustering y visualizar:

set.seed(123)
kmeans_result <- kmeans(data_scaled, centers = 6, nstart = 25)

# Visualizar clústeres en el espacio PCA
fviz_cluster(kmeans_result, data = data_scaled,
             ellipse.type = "norm",
             geom = "point",
             palette = "jco",
             ggtheme = theme_minimal()) +
  theme_nyc()

Vemos que no se obtiene nada claro ni interesante para analizar. Vamos a probar con otro tipo de clustering, el clustering jerárquico.

4.2.1.4 Clustering jerarquico

# Calcular la matriz de distancias
dist_matrix <- dist(data_scaled, method = "euclidean")

# Realizar el clustering jerárquico
hc <- hclust(dist_matrix, method = "complete") # --> Mirar a ver que método es mejor

# Visualizar el dendrograma
plot(hc, labels = FALSE, hang = -1, main = "Dendrograma de Clustering Jerárquico")

# Cortamos:
rect.hclust(hc, k = 3, border = "red")

clusters <- cutree(hc, k = 3)
table(clusters)

## clusters
##     1     2     3 
## 19988     9     3

Dado que tal y como hemos calculado el clúster jerárquico no obtenemos un resultado interesante, vamos a probar la alternativa de agrupar los datos por barrios.

4.2.1.5 Agrupar los datos por barrios y clusterizar:

# Sólo con el total de muertos y heridos:
data_borough_sum <- data_sampled |> 
  group_by(BOROUGH) |> 
  summarise(
    total_injured = sum(NUM_PERSONS_INJURED, na.rm = TRUE),
    total_killed = sum(NUM_PERSONS_KILLED, na.rm = TRUE)
  ) |> 
  na.omit()

Escalar los datos:

data_scaled <- data_borough_sum |> 
  dplyr::select(-BOROUGH) |> 
  scale()

Clustering jerárquico:

dist_matrix <- dist(data_scaled)
hc <- hclust(dist_matrix, method = "complete")

# Dendrograma
dend <- as.dendrogram(hc)
dend <- color_branches(dend, k = 3)
dend <- set(dend, "labels", data_borough_sum$BOROUGH[hc$order])

factoextra::fviz_dend(
  dend,
  k = 3,
  horiz = TRUE,
  rect = TRUE,
  rect_border = "#1F3B73",
  rect_fill = FALSE,
  main = "Dendrograma de los Distritos",
  cex = 0.7,
  color_labels_by_k = TRUE,
) + 
theme_nyc()

El análisis de clúster jerárquico aplicado a los distritos de Nueva York ha dado como resultado una agrupación en tres clústeres bien diferenciados: por un lado, se agrupan el Bronx y Manhattan; por otro, Brooklyn y Queens; y finalmente, Staten Island queda separado en un grupo propio. Esta clasificación cobra sentido cuando se analiza en función de la densidad de población de cada distrito.

Manhattan, con 27.330 habitantes por kilómetro cuadrado, y el Bronx, con 12.831, presentan las mayores densidades de población entre los distritos analizados. Esta alta concentración de personas suele asociarse a una mayor exposición al tráfico, lo que incrementa la probabilidad de accidentes y, por tanto, de víctimas. Esto justificaría que ambos formen parte del mismo grupo, ya que comparten patrones similares tanto en cuanto a intensidad del tráfico como al número total de heridos y fallecidos.

En un segundo grupo aparecen Brooklyn y Queens, con densidades de población intermedias (14.400 y 8.104,7 habitantes por km² respectivamente). Aunque Brooklyn supera al Bronx en densidad, su comportamiento en términos de accidentes se asemeja más al de Queens, probablemente por la estructura urbana de ambos: grandes áreas residenciales, presencia de autopistas y una combinación de zonas densas y otras más dispersas. Esta mezcla da lugar a un número total de víctimas por accidentes más moderado, lo que explica su asociación en un mismo clúster.

Por último, Staten Island se mantiene como un grupo aislado, con una densidad de apenas 1.782,2 habitantes por km². Este distrito, claramente menos poblado y urbanizado, muestra un patrón de accidentes significativamente distinto, con un número mucho menor de víctimas en comparación con el resto. Su baja densidad y características geográficas más suburbanas justifican su separación en un clúster independiente.

Esta clasificación por clústeres, por tanto, se alinea coherentemente con la densidad de población de cada distrito y permite explicar diferencias en el comportamiento de los accidentes de tráfico a lo largo de la ciudad.

Vamos a añadir un gráfico de barras comparando las variables entre los clústeres:

1. Asignamos los clúster al dataframe:

# Cortar el dendrograma para obtener 3 clústeres
clusters <- cutree(hc, k = 3)

# Añadir la columna de clúster al dataframe original
data_borough_sum$cluster <- as.factor(clusters)

2. Convertimos a formato largo para ggplot:

data_borough_long <- data_borough_sum |> 
  pivot_longer(cols = c(total_injured, total_killed), 
               names_to = "Variable", 
               values_to = "Total")

3. Creamos el gráfico de barras:

ggplot(data_borough_long, aes(x = Variable, y = Total, fill = cluster)) +
  geom_col(position = "dodge") +
  facet_wrap(~BOROUGH) +
  labs(title = "Comparación de heridos y muertos por clúster y barrio",
       x = "Variable", y = "Total", fill = "Clúster") +
  theme_nyc()

Este gráfico demuestra que nuestra hipótesis anterior tiene sentido, ya que el clúster es capaz de agrupar los distintos distritos de una manera muy correcta en función del número total de heridos y muertos.

Por último, añadiremos el shilouette para comprobar si nuestra disposición del clúster es buena o no:

clusters <- cutree(hc, k = 3)

dist_matrix <- dist(data_borough_sum[, c("total_injured", "total_killed")])

sil <- silhouette(clusters, dist_matrix)
rownames(sil) <- data_borough_sum$BOROUGH

fviz_silhouette(sil) + 
  ggtitle(paste("Índice de Silueta para k =", k)) +
  theme_nyc()

##   cluster size ave.sil.width
## 1       1    2          0.83
## 2       2    2          0.91
## 3       3    1          0.00

El análisis de la silueta respalda en buena medida la clasificación obtenida. Los grupos formados por el Bronx junto a Manhattan, y Brooklyn junto a Queens, muestran una estructura muy sólida: sus valores de silueta son claramente altos, lo que indica que los distritos dentro de cada clúster se parecen mucho entre sí y están bien separados de los demás. Esto refuerza la idea de que ambos grupos tienen sentido, tanto desde un punto de vista numérico como contextual.

Sin embargo, Staten Island, que ha quedado como un clúster independiente, presenta una silueta bastante más baja, de 0. Esto quiere decir que, aunque por sus características es razonable considerarlo un caso aparte, los datos no lo separan de forma tan clara como ocurre con los otros clústeres. Aun así, teniendo en cuenta su baja densidad de población y su configuración más suburbana, parece justificado mantenerlo como un grupo propio. Este contraste entre lo que dicen los números y lo que observamos en la realidad muestra lo útil que es combinar el análisis estadístico con la interpretación del contexto para entender mejor los resultados.

4.2.2 Análisis Cluster Jerárquico variable ‘CAUSE’:

En cuanto a la variable CAUSE, correspondiente a los factores contribuyentes al accidente, nos ha parecido interesante aplicar también clusterización ya que tenemos muchas de ellas y de esta forma veremos si conseguimos reducirlas a un número razonable. Para ello, vamos a usar el algoritmo K-means, que nos permite agrupar los diferentes factores en diferentes grupos.

Establecemos similitudes entre los nombres de las causas, mediante una matriz de distancias que representa la distancia entre las diferentes causas según su similitud. Para ello lo hacemos mediante el metodo de “cosine” ya que se basa en frecuencia de caracteres y es bueno si las cadenas comparten subcadenas similares.

# Crear matriz de distancias entre nombres de causas
dist_matrix <- stringdistmatrix(causes_cleaned$CAUSE, causes_cleaned$CAUSE, method = "cosine")  

rownames(dist_matrix) <- causes_cleaned$CAUSE
colnames(dist_matrix) <- causes_cleaned$CAUSE

Para aplicar clustering jerarquico, tomamos como método el “complete”, que es una medida más conservadora. También hacemos un test de gap para que nos proporcione el mejor valor de k.

hc <- hclust(as.dist(dist_matrix), method = "complete")

fviz_nbclust(as.matrix(dist_matrix), FUN = hcut, method = "gap_stat", diss = TRUE)

El test nos dice que el valor óptimo es k = 2. Ahora vamos a realizar el clustering jerárquico y graficar los resultados, junto a la silueta para ver que tal realiza la clusterización:

k <- 2
fviz_dend(x = hc,
          k = k,
          k_colors = c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", 
                       "#9467bd", "#8c564b",  "#e377c2", "#7f7f7f"),
          color_labels_by_k = TRUE,
          rect = TRUE,
          rect_fill = TRUE,
          cex = 0.5,
          main = "Dendrograma - ward",
          xlab = "observaciones",
          ylab = "distancia",
          sub = "") +
  theme_nyc()

clusters <- cutree(hc, k = k)

# Calcular silueta
sil <- silhouette(clusters, dist(as.dist(dist_matrix)))

# Graficar con estilo personalizado
fviz_silhouette(sil) +
  ggtitle(paste("Índice de Silueta para k =", k)) +
  theme_nyc()

##   cluster size ave.sil.width
## 1       1   31          0.11
## 2       2   21          0.16

Aunque el método del Gap statistic sugiere que k = 2 es óptimo, la silueta indica que la calidad del agrupamiento es baja. Esto sugiere que las causas no están bien separadas en dos grupos distintos. Sin embargo, esto también puede ser un reflejo de la complejidad y diversidad de las causas de los accidentes, lo que hace difícil encontrar una separación clara entre ellas.

Ahora vamos a probar con un valor de k = 9, que en el test de gap nos daba el valor más alto, a ver si conseguimos una mejor separación entre las causas. Esto nos permitirá ver si hay patrones más claros en los datos y si podemos agrupar las causas de manera más efectiva.

k <- 9
fviz_dend(x = hc,
          k = k,
          k_colors = c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b",  "#e377c2", "#7f7f7f"),
          color_labels_by_k = TRUE,
          rect = TRUE,
          rect_fill = TRUE,
          cex = 0.5,
          main = "Dendrograma - ward",
          xlab = "observaciones",
          ylab = "distancia",
          sub = "") +
  theme_nyc()

clusters <- cutree(hc, k = k)

# Calcular silueta
sil <- silhouette(clusters, dist(as.dist(dist_matrix)))

# Graficar con estilo personalizado
fviz_silhouette(sil) +
  ggtitle(paste("Índice de Silueta para k =", k)) +
  theme_nyc()

##   cluster size ave.sil.width
## 1       1   19          0.04
## 2       2    3          0.28
## 3       3    4          0.11
## 4       4    8         -0.01
## 5       5    4          0.10
## 6       6    9          0.24
## 7       7    2          0.03
## 8       8    1          0.00
## 9       9    2          0.33

Para un valor de k = 9 se pueden observar cosas interesantes. A pesar de que haya valores negativos, hay algunos cluster bastante bien agrupados con causas bastante similares, como por ejemplo la última rama del endograma, ya que reune todas las causas relacionades con defectos, ya sea en el pavimento, en los frenos…

4.3 Análisis discriminante

Vamos a predecir el distrito al que pertenece un accidente. Para ello vamos a aseguir los siguientes pasos:

4.3.1 1. Seleccionar las variables predictoras y agrupar categorías poco frecuentes:

# Agrupar las 10 categorías más frecuentes de VEHICLE_1
top_veh_types <- data_sampled |> 
  count(VEHICLE_1, sort = TRUE) |> 
  slice_head(n = 10) |> 
  pull(VEHICLE_1)

data_sampled <- data_sampled |> 
  mutate(VEHICLE_1 = ifelse(VEHICLE_1 %in% top_veh_types, VEHICLE_1, "OTROS"))

# Agrupar las 10 categorías más frecuentes de VEHICLE_2
top_veh_types <- data_sampled |> 
  count(VEHICLE_2, sort = TRUE) |> 
  slice_head(n = 10) |> 
  pull(VEHICLE_2)

data_sampled <- data_sampled |> 
  mutate(VEHICLE_2 = ifelse(VEHICLE_2 %in% top_veh_types, VEHICLE_2, "OTROS"))

# Agrupar las 10 categorías más frecuentes de CAUSE
top_causes <- data_sampled |> 
  count(CAUSE, sort = TRUE) |> 
  slice_head(n = 10) |> 
  pull(CAUSE)

data_sampled <- data_sampled |> 
  mutate(CAUSE = ifelse(CAUSE %in% top_causes, CAUSE, "OTROS"))

Seleccionamos las variables:

data_lda <- data_sampled |> 
  dplyr::select(BOROUGH,
         NUM_PERSONS_INJURED, NUM_PERSONS_KILLED, LONGITUDE, LATITUDE,
         NUM_PEDESTRIANS_INJURED, NUM_PEDESTRIANS_KILLED,
         NUM_CYCLIST_INJURED, NUM_CYCLIST_KILLED,
         NUM_MOTORIST_INJURED, NUM_MOTORIST_KILLED,
         VEHICLE_1, VEHICLE_2, CAUSE) |>
  na.omit()

4.3.2 2. Convertir las variables categóricas en dummies:

data_lda_dummy <- dummyVars(BOROUGH ~ ., data = data_lda)
data_lda_transformed <- predict(data_lda_dummy, newdata = data_lda)
data_lda_final <- data.frame(BOROUGH = data_lda$BOROUGH, data_lda_transformed)

4.3.3 3. Crear un conjunto de entrenamiento y test:

set.seed(123)
train_index <- createDataPartition(data_lda_final$BOROUGH, p = 0.7, list = FALSE)
train_data <- data_lda_final[train_index, ]
test_data <- data_lda_final[-train_index, ]

4.3.4 4. Aplicamos el análisis discriminante:

# Aplicar LDA
lda_model <- lda(BOROUGH ~ ., data = train_data)

## Warning in lda.default(x, grouping, ...): variables are collinear

lda_model

## Call:
## lda(BOROUGH ~ ., data = train_data)
## 
## Prior probabilities of groups:
##         BRONX      BROOKLYN     MANHATTAN        QUEENS STATEN ISLAND 
##    0.13733752    0.29502928    0.24446508    0.28231681    0.04085131 
## 
## Group means:
##               NUM_PERSONS_INJURED NUM_PERSONS_KILLED LONGITUDE LATITUDE
## BRONX                   0.3510140       0.0005200208 -73.88458 40.84716
## BROOKLYN                0.3464052       0.0019365771 -73.95252 40.65879
## MANHATTAN               0.2290389       0.0002921414 -73.97603 40.76271
## QUEENS                  0.3518846       0.0012648621 -73.83243 40.71934
## STATEN ISLAND           0.3986014       0.0017482517 -74.12373 40.58424
##               NUM_PEDESTRIANS_INJURED NUM_PEDESTRIANS_KILLED
## BRONX                      0.07332293           0.0000000000
## BROOKLYN                   0.06463326           0.0007262164
## MANHATTAN                  0.06777680           0.0002921414
## QUEENS                     0.05008854           0.0005059449
## STATEN ISLAND              0.04895105           0.0000000000
##               NUM_CYCLIST_INJURED NUM_CYCLIST_KILLED NUM_MOTORIST_INJURED
## BRONX                  0.02600104       0.0000000000            0.2516901
## BROOKLYN               0.04720407       0.0004841443            0.2345679
## MANHATTAN              0.03739410       0.0000000000            0.1238680
## QUEENS                 0.02833291       0.0000000000            0.2734632
## STATEN ISLAND          0.01048951       0.0000000000            0.3391608
##               NUM_MOTORIST_KILLED VEHICLE_14.DR.SEDAN VEHICLE_1BOX.TRUCK
## BRONX                0.0005200208          0.01716069        0.010920437
## BROOKLYN             0.0007262164          0.01767127        0.017671266
## MANHATTAN            0.0000000000          0.01635992        0.024247736
## QUEENS               0.0007589173          0.01947888        0.009359980
## STATEN ISLAND        0.0017482517          0.01573427        0.001748252
##               VEHICLE_1BUS VEHICLE_1OTROS VEHICLE_1PASSENGER.VEHICLE
## BRONX           0.01976079     0.08736349                  0.1378055
## BROOKLYN        0.01815541     0.08569354                  0.1508109
## MANHATTAN       0.02453988     0.12094654                  0.1560035
## QUEENS          0.01062484     0.05919555                  0.1573488
## STATEN ISLAND   0.01398601     0.03321678                  0.1538462
##               VEHICLE_1PICK.UP.TRUCK VEHICLE_1SEDAN
## BRONX                     0.02080083      0.3416537
## BROOKLYN                  0.02227064      0.3367223
## MANHATTAN                 0.02395559      0.2196903
## QUEENS                    0.02099671      0.3303820
## STATEN ISLAND             0.03496503      0.4160839
##               VEHICLE_1SPORT.UTILITY...STATION.WAGON
## BRONX                                     0.05824233
## BROOKLYN                                  0.06197047
## MANHATTAN                                 0.07420391
## QUEENS                                    0.07715659
## STATEN ISLAND                             0.08041958
##               VEHICLE_1STATION.WAGON.SPORT.UTILITY.VEHICLE VEHICLE_1TAXI
## BRONX                                            0.2652106   0.030681227
## BROOKLYN                                         0.2578068   0.016218833
## MANHATTAN                                        0.1735320   0.137014315
## QUEENS                                           0.2914242   0.016190235
## STATEN ISLAND                                    0.2342657   0.005244755
##               VEHICLE_1VAN VEHICLE_2 VEHICLE_2BIKE VEHICLE_2BUS VEHICLE_2OTROS
## BRONX          0.010400416 0.2033281    0.01456058   0.01300052     0.11232449
## BROOKLYN       0.015008473 0.1745340    0.02832244   0.01597676     0.12611958
## MANHATTAN      0.029506281 0.1504528    0.02453988   0.02337131     0.16798130
## QUEENS         0.007842145 0.1611434    0.01745510   0.01467240     0.08651657
## STATEN ISLAND  0.010489510 0.2097902    0.01048951   0.01223776     0.05594406
##               VEHICLE_2PASSENGER.VEHICLE VEHICLE_2PICK.UP.TRUCK VEHICLE_2SEDAN
## BRONX                          0.1050442             0.02132085      0.2792512
## BROOKLYN                       0.1099008             0.01791334      0.2394093
## MANHATTAN                      0.1194858             0.02600058      0.1688577
## QUEENS                         0.1300278             0.02302049      0.2544903
## STATEN ISLAND                  0.1293706             0.02797203      0.2709790
##               VEHICLE_2SPORT.UTILITY...STATION.WAGON
## BRONX                                     0.04524181
## BROOKLYN                                  0.05422416
## MANHATTAN                                 0.06485539
## QUEENS                                    0.06121933
## STATEN ISLAND                             0.06818182
##               VEHICLE_2STATION.WAGON.SPORT.UTILITY.VEHICLE VEHICLE_2TAXI
## BRONX                                            0.1684867   0.013520541
## BROOKLYN                                         0.2023723   0.012103607
## MANHATTAN                                        0.1142273   0.117732983
## QUEENS                                           0.2203390   0.011383759
## STATEN ISLAND                                    0.2010490   0.001748252
##               VEHICLE_2UNKNOWN CAUSEBACKING.UNSAFELY
## BRONX               0.02392096            0.07384295
## BROOKLYN            0.01912370            0.07237957
## MANHATTAN           0.02249489            0.05054046
## QUEENS              0.01973185            0.07437389
## STATEN ISLAND       0.01223776            0.06293706
##               CAUSEDRIVER.INATTENTION.DISTRACTION
## BRONX                                   0.2693708
## BROOKLYN                                0.2931494
## MANHATTAN                               0.3324569
## QUEENS                                  0.3313939
## STATEN ISLAND                           0.3583916
##               CAUSEFAILURE.TO.YIELD.RIGHT.OF.WAY CAUSEFATIGUED.DROWSY
## BRONX                                 0.08476339           0.01508060
## BROOKLYN                              0.10990075           0.03243767
## MANHATTAN                             0.06982179           0.04060765
## QUEENS                                0.13028080           0.02403238
## STATEN ISLAND                         0.12587413           0.02797203
##               CAUSEFOLLOWING.TOO.CLOSELY CAUSEOTHER.VEHICULAR CAUSEOTROS
## BRONX                         0.04732189           0.08632345  0.2542902
## BROOKLYN                      0.05664488           0.03582668  0.2444929
## MANHATTAN                     0.04557406           0.08735028  0.2255332
## QUEENS                        0.05211232           0.02428535  0.2135087
## STATEN ISLAND                 0.04720280           0.02972028  0.2587413
##               CAUSEPASSING.OR.LANE.USAGE.IMPROPER CAUSEPASSING.TOO.CLOSELY
## BRONX                                  0.03224129               0.04784191
## BROOKLYN                               0.04139434               0.05204551
## MANHATTAN                              0.04206836               0.03856266
## QUEENS                                 0.04578801               0.02959777
## STATEN ISLAND                          0.02272727               0.01748252
##               CAUSETRAFFIC.CONTROL.DISREGARDED CAUSETURNING.IMPROPERLY
## BRONX                               0.03744150              0.05148206
## BROOKLYN                            0.03485839              0.02687001
## MANHATTAN                           0.01694420              0.05054046
## QUEENS                              0.04123451              0.03339236
## STATEN ISLAND                       0.02272727              0.02622378
## 
## Coefficients of linear discriminants:
##                                                        LD1           LD2
## NUM_PERSONS_INJURED                           1.787585e-02  2.459970e-03
## NUM_PERSONS_KILLED                            8.958580e-02 -1.895147e-01
## LONGITUDE                                    -1.651256e+01  1.857395e+01
## LATITUDE                                     -1.920146e+01 -1.859816e+01
## NUM_PEDESTRIANS_INJURED                      -8.025201e-02 -1.183892e-01
## NUM_PEDESTRIANS_KILLED                        5.961387e-01  1.456901e-01
## NUM_CYCLIST_INJURED                           8.266605e-02  2.497290e-01
## NUM_CYCLIST_KILLED                            1.131411e+00 -3.037183e-01
## NUM_MOTORIST_INJURED                          2.560329e-02  5.337511e-04
## NUM_MOTORIST_KILLED                          -5.511250e-01 -4.121097e-01
## VEHICLE_14.DR.SEDAN                          -1.075601e-03  9.498872e-03
## VEHICLE_1BOX.TRUCK                           -1.010474e-04  3.383678e-02
## VEHICLE_1BUS                                  2.915018e-02 -2.201628e-01
## VEHICLE_1OTROS                               -9.278053e-02  2.554645e-02
## VEHICLE_1PASSENGER.VEHICLE                    3.457094e-02  3.458676e-02
## VEHICLE_1PICK.UP.TRUCK                        6.993791e-02 -1.258601e-02
## VEHICLE_1SEDAN                                6.277424e-02 -4.557795e-03
## VEHICLE_1SPORT.UTILITY...STATION.WAGON       -4.918220e-02  7.195133e-02
## VEHICLE_1STATION.WAGON.SPORT.UTILITY.VEHICLE -4.095841e-03  8.968837e-02
## VEHICLE_1TAXI                                -1.971117e-01 -4.661447e-01
## VEHICLE_1VAN                                 -1.948667e-02 -2.229110e-01
## VEHICLE_2                                     3.899696e-02  3.262172e-02
## VEHICLE_2BIKE                                -1.115658e-01 -1.992471e-02
## VEHICLE_2BUS                                 -1.058006e-02 -1.209988e-01
## VEHICLE_2OTROS                               -1.047567e-01 -3.411731e-02
## VEHICLE_2PASSENGER.VEHICLE                    3.091809e-02  9.407475e-04
## VEHICLE_2PICK.UP.TRUCK                        3.155406e-02  4.080690e-02
## VEHICLE_2SEDAN                                3.179511e-02 -1.081030e-02
## VEHICLE_2SPORT.UTILITY...STATION.WAGON       -3.958150e-02 -1.592596e-02
## VEHICLE_2STATION.WAGON.SPORT.UTILITY.VEHICLE  5.063315e-02  1.145683e-01
## VEHICLE_2TAXI                                -1.564345e-01 -4.243846e-01
## VEHICLE_2UNKNOWN                             -7.232373e-02  4.168859e-02
## CAUSEBACKING.UNSAFELY                         7.890740e-03  9.330540e-02
## CAUSEDRIVER.INATTENTION.DISTRACTION          -5.343581e-02 -2.607940e-03
## CAUSEFAILURE.TO.YIELD.RIGHT.OF.WAY            5.025285e-02  1.247055e-01
## CAUSEFATIGUED.DROWSY                          6.051003e-03  3.079434e-02
## CAUSEFOLLOWING.TOO.CLOSELY                    2.971536e-03  1.321781e-01
## CAUSEOTHER.VEHICULAR                          4.752801e-02 -4.499541e-01
## CAUSEOTROS                                    1.915428e-02 -2.039519e-02
## CAUSEPASSING.OR.LANE.USAGE.IMPROPER           1.437565e-02  8.292922e-02
## CAUSEPASSING.TOO.CLOSELY                      7.461936e-02  4.373198e-02
## CAUSETRAFFIC.CONTROL.DISREGARDED             -3.884305e-03  2.328326e-02
## CAUSETURNING.IMPROPERLY                      -8.318955e-02 -1.177570e-01
##                                                       LD3         LD4
## NUM_PERSONS_INJURED                           0.126569460 -0.04897919
## NUM_PERSONS_KILLED                            0.356742104 -0.60594932
## LONGITUDE                                    -3.629841563 -0.57283358
## LATITUDE                                      3.629398298  0.46886749
## NUM_PEDESTRIANS_INJURED                      -0.486045425 -1.28693880
## NUM_PEDESTRIANS_KILLED                       -1.545600199 -2.45486542
## NUM_CYCLIST_INJURED                           0.125924546 -1.41680835
## NUM_CYCLIST_KILLED                            1.467479994 -8.82469948
## NUM_MOTORIST_INJURED                          0.205699850  0.25034543
## NUM_MOTORIST_KILLED                           1.505497380  2.83751975
## VEHICLE_14.DR.SEDAN                           0.345317095  0.49141930
## VEHICLE_1BOX.TRUCK                           -1.128507000 -1.41132376
## VEHICLE_1BUS                                 -0.247052371 -0.92650842
## VEHICLE_1OTROS                               -0.577474720 -0.88761994
## VEHICLE_1PASSENGER.VEHICLE                    0.072867293  0.02435347
## VEHICLE_1PICK.UP.TRUCK                        0.066233576  0.56956564
## VEHICLE_1SEDAN                                0.486013887  0.06167821
## VEHICLE_1SPORT.UTILITY...STATION.WAGON        0.008810391  0.64933911
## VEHICLE_1STATION.WAGON.SPORT.UTILITY.VEHICLE  0.335574860  0.01886118
## VEHICLE_1TAXI                                -2.315320956  0.79462598
## VEHICLE_1VAN                                 -1.152607498 -0.69396991
## VEHICLE_2                                     0.504261559  0.29596475
## VEHICLE_2BIKE                                -0.754308798 -0.69340399
## VEHICLE_2BUS                                 -0.692906603 -0.14332392
## VEHICLE_2OTROS                               -0.567299320 -0.99253160
## VEHICLE_2PASSENGER.VEHICLE                    0.133600124  0.35740797
## VEHICLE_2PICK.UP.TRUCK                       -0.176402601  0.91210954
## VEHICLE_2SEDAN                                0.372468752 -0.06517169
## VEHICLE_2SPORT.UTILITY...STATION.WAGON       -0.055574083  0.39265543
## VEHICLE_2STATION.WAGON.SPORT.UTILITY.VEHICLE  0.260473050 -0.03865669
## VEHICLE_2TAXI                                -2.882418594  0.78131784
## VEHICLE_2UNKNOWN                              0.130305172 -0.41714658
## CAUSEBACKING.UNSAFELY                         0.265105080 -0.23925456
## CAUSEDRIVER.INATTENTION.DISTRACTION          -0.196050835  0.61991069
## CAUSEFAILURE.TO.YIELD.RIGHT.OF.WAY            0.204475208  0.64805168
## CAUSEFATIGUED.DROWSY                         -0.417006164 -0.63778323
## CAUSEFOLLOWING.TOO.CLOSELY                   -0.122753060 -0.40621896
## CAUSEOTHER.VEHICULAR                         -0.158396839 -0.54135516
## CAUSEOTROS                                    0.206183893 -0.19004590
## CAUSEPASSING.OR.LANE.USAGE.IMPROPER          -0.464917485 -0.10482899
## CAUSEPASSING.TOO.CLOSELY                     -0.042475087 -1.92263761
## CAUSETRAFFIC.CONTROL.DISREGARDED              0.419404979 -0.32674493
## CAUSETURNING.IMPROPERLY                       0.058206646  0.25274079
## 
## Proportion of trace:
##    LD1    LD2    LD3    LD4 
## 0.6349 0.3492 0.0137 0.0021

4.3.5 5. Evaluamos el modelo:

predictions <- predict(lda_model, newdata = test_data)

# Matriz de confusión
table(Predicho = predictions$class, Real = test_data$BOROUGH)

##                Real
## Predicho        BRONX BROOKLYN MANHATTAN QUEENS STATEN ISLAND
##   BRONX           801        0       107      1             0
##   BROOKLYN          0     1643        36     45             2
##   MANHATTAN        21       83      1324    202             0
##   QUEENS            1       42         0   1446             0
##   STATEN ISLAND     0        2         0      0           242

4.3.6 6. Visualizar las dos primeras funciones discriminantes:

lda_df <- data.frame(predictions$x, BOROUGH = test_data$BOROUGH)

ggplot(lda_df, aes(x = LD1, y = LD2, color = BOROUGH)) +
  geom_point(alpha = 0.5) +
  labs(title = "Proyección LDA: Accidentes por distrito",
       x = "Función discriminante 1",
       y = "Función discriminante 2") +
  theme_nyc()

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

Se observa que el modelo obtiene buenos resultados únicamente cuando se incluyen las variables LONGITUDE y LATITUDE. Sin embargo, esto introduce una especie de “trampa”, ya que esas variables están muy relacionadas con la ubicación exacta del accidente, lo que facilita en exceso la predicción. Por ello, para evaluar realmente la capacidad del modelo, resulta más apropiado repetir el análisis excluyendo estas variables y observar cómo se comporta con información menos directa.

Vamos a comprobar entonces qué pasaría si realizaramos el mismo análisis sin esas variables:

4.3.7 1.2 Selección de variables:

Seleccionamos las variables, excepto LATITUDE y LONGITUDE:

data_lda <- data_sampled |> 
  dplyr::select(BOROUGH,
         NUM_PERSONS_INJURED, NUM_PERSONS_KILLED,
         NUM_PEDESTRIANS_INJURED, NUM_PEDESTRIANS_KILLED,
         NUM_CYCLIST_INJURED, NUM_CYCLIST_KILLED,
         NUM_MOTORIST_INJURED, NUM_MOTORIST_KILLED,
         VEHICLE_1, VEHICLE_2, CAUSE) |>
  na.omit()

4.3.8 2.2 Convertir las variables categóricas en dummies:

data_lda_dummy <- dummyVars(BOROUGH ~ ., data = data_lda)
data_lda_transformed <- predict(data_lda_dummy, newdata = data_lda)
data_lda_final <- data.frame(BOROUGH = data_lda$BOROUGH, data_lda_transformed)

4.3.9 3.2 Crear un conjunto de entrenamiento y test:

set.seed(123)
train_index <- createDataPartition(data_lda_final$BOROUGH, p = 0.7, list = FALSE)
train_data <- data_lda_final[train_index, ]
test_data <- data_lda_final[-train_index, ]

4.3.10 4.2 Aplicamos el análisis discriminante:

# Aplicar LDA
lda_model <- lda(BOROUGH ~ ., data = train_data)

## Warning in lda.default(x, grouping, ...): variables are collinear

lda_model

## Call:
## lda(BOROUGH ~ ., data = train_data)
## 
## Prior probabilities of groups:
##         BRONX      BROOKLYN     MANHATTAN        QUEENS STATEN ISLAND 
##    0.13733752    0.29502928    0.24446508    0.28231681    0.04085131 
## 
## Group means:
##               NUM_PERSONS_INJURED NUM_PERSONS_KILLED NUM_PEDESTRIANS_INJURED
## BRONX                   0.3510140       0.0005200208              0.07332293
## BROOKLYN                0.3464052       0.0019365771              0.06463326
## MANHATTAN               0.2290389       0.0002921414              0.06777680
## QUEENS                  0.3518846       0.0012648621              0.05008854
## STATEN ISLAND           0.3986014       0.0017482517              0.04895105
##               NUM_PEDESTRIANS_KILLED NUM_CYCLIST_INJURED NUM_CYCLIST_KILLED
## BRONX                   0.0000000000          0.02600104       0.0000000000
## BROOKLYN                0.0007262164          0.04720407       0.0004841443
## MANHATTAN               0.0002921414          0.03739410       0.0000000000
## QUEENS                  0.0005059449          0.02833291       0.0000000000
## STATEN ISLAND           0.0000000000          0.01048951       0.0000000000
##               NUM_MOTORIST_INJURED NUM_MOTORIST_KILLED VEHICLE_14.DR.SEDAN
## BRONX                    0.2516901        0.0005200208          0.01716069
## BROOKLYN                 0.2345679        0.0007262164          0.01767127
## MANHATTAN                0.1238680        0.0000000000          0.01635992
## QUEENS                   0.2734632        0.0007589173          0.01947888
## STATEN ISLAND            0.3391608        0.0017482517          0.01573427
##               VEHICLE_1BOX.TRUCK VEHICLE_1BUS VEHICLE_1OTROS
## BRONX                0.010920437   0.01976079     0.08736349
## BROOKLYN             0.017671266   0.01815541     0.08569354
## MANHATTAN            0.024247736   0.02453988     0.12094654
## QUEENS               0.009359980   0.01062484     0.05919555
## STATEN ISLAND        0.001748252   0.01398601     0.03321678
##               VEHICLE_1PASSENGER.VEHICLE VEHICLE_1PICK.UP.TRUCK VEHICLE_1SEDAN
## BRONX                          0.1378055             0.02080083      0.3416537
## BROOKLYN                       0.1508109             0.02227064      0.3367223
## MANHATTAN                      0.1560035             0.02395559      0.2196903
## QUEENS                         0.1573488             0.02099671      0.3303820
## STATEN ISLAND                  0.1538462             0.03496503      0.4160839
##               VEHICLE_1SPORT.UTILITY...STATION.WAGON
## BRONX                                     0.05824233
## BROOKLYN                                  0.06197047
## MANHATTAN                                 0.07420391
## QUEENS                                    0.07715659
## STATEN ISLAND                             0.08041958
##               VEHICLE_1STATION.WAGON.SPORT.UTILITY.VEHICLE VEHICLE_1TAXI
## BRONX                                            0.2652106   0.030681227
## BROOKLYN                                         0.2578068   0.016218833
## MANHATTAN                                        0.1735320   0.137014315
## QUEENS                                           0.2914242   0.016190235
## STATEN ISLAND                                    0.2342657   0.005244755
##               VEHICLE_1VAN VEHICLE_2 VEHICLE_2BIKE VEHICLE_2BUS VEHICLE_2OTROS
## BRONX          0.010400416 0.2033281    0.01456058   0.01300052     0.11232449
## BROOKLYN       0.015008473 0.1745340    0.02832244   0.01597676     0.12611958
## MANHATTAN      0.029506281 0.1504528    0.02453988   0.02337131     0.16798130
## QUEENS         0.007842145 0.1611434    0.01745510   0.01467240     0.08651657
## STATEN ISLAND  0.010489510 0.2097902    0.01048951   0.01223776     0.05594406
##               VEHICLE_2PASSENGER.VEHICLE VEHICLE_2PICK.UP.TRUCK VEHICLE_2SEDAN
## BRONX                          0.1050442             0.02132085      0.2792512
## BROOKLYN                       0.1099008             0.01791334      0.2394093
## MANHATTAN                      0.1194858             0.02600058      0.1688577
## QUEENS                         0.1300278             0.02302049      0.2544903
## STATEN ISLAND                  0.1293706             0.02797203      0.2709790
##               VEHICLE_2SPORT.UTILITY...STATION.WAGON
## BRONX                                     0.04524181
## BROOKLYN                                  0.05422416
## MANHATTAN                                 0.06485539
## QUEENS                                    0.06121933
## STATEN ISLAND                             0.06818182
##               VEHICLE_2STATION.WAGON.SPORT.UTILITY.VEHICLE VEHICLE_2TAXI
## BRONX                                            0.1684867   0.013520541
## BROOKLYN                                         0.2023723   0.012103607
## MANHATTAN                                        0.1142273   0.117732983
## QUEENS                                           0.2203390   0.011383759
## STATEN ISLAND                                    0.2010490   0.001748252
##               VEHICLE_2UNKNOWN CAUSEBACKING.UNSAFELY
## BRONX               0.02392096            0.07384295
## BROOKLYN            0.01912370            0.07237957
## MANHATTAN           0.02249489            0.05054046
## QUEENS              0.01973185            0.07437389
## STATEN ISLAND       0.01223776            0.06293706
##               CAUSEDRIVER.INATTENTION.DISTRACTION
## BRONX                                   0.2693708
## BROOKLYN                                0.2931494
## MANHATTAN                               0.3324569
## QUEENS                                  0.3313939
## STATEN ISLAND                           0.3583916
##               CAUSEFAILURE.TO.YIELD.RIGHT.OF.WAY CAUSEFATIGUED.DROWSY
## BRONX                                 0.08476339           0.01508060
## BROOKLYN                              0.10990075           0.03243767
## MANHATTAN                             0.06982179           0.04060765
## QUEENS                                0.13028080           0.02403238
## STATEN ISLAND                         0.12587413           0.02797203
##               CAUSEFOLLOWING.TOO.CLOSELY CAUSEOTHER.VEHICULAR CAUSEOTROS
## BRONX                         0.04732189           0.08632345  0.2542902
## BROOKLYN                      0.05664488           0.03582668  0.2444929
## MANHATTAN                     0.04557406           0.08735028  0.2255332
## QUEENS                        0.05211232           0.02428535  0.2135087
## STATEN ISLAND                 0.04720280           0.02972028  0.2587413
##               CAUSEPASSING.OR.LANE.USAGE.IMPROPER CAUSEPASSING.TOO.CLOSELY
## BRONX                                  0.03224129               0.04784191
## BROOKLYN                               0.04139434               0.05204551
## MANHATTAN                              0.04206836               0.03856266
## QUEENS                                 0.04578801               0.02959777
## STATEN ISLAND                          0.02272727               0.01748252
##               CAUSETRAFFIC.CONTROL.DISREGARDED CAUSETURNING.IMPROPERLY
## BRONX                               0.03744150              0.05148206
## BROOKLYN                            0.03485839              0.02687001
## MANHATTAN                           0.01694420              0.05054046
## QUEENS                              0.04123451              0.03339236
## STATEN ISLAND                       0.02272727              0.02622378
## 
## Coefficients of linear discriminants:
##                                                       LD1         LD2
## NUM_PERSONS_INJURED                          -0.109402051  0.07857783
## NUM_PERSONS_KILLED                           -0.568437652 -0.38274147
## NUM_PEDESTRIANS_INJURED                       0.589091789  0.87269788
## NUM_PEDESTRIANS_KILLED                        0.542395413 -1.03285596
## NUM_CYCLIST_INJURED                          -0.400485523  0.22072651
## NUM_CYCLIST_KILLED                           -1.644881635  3.58914279
## NUM_MOTORIST_INJURED                         -0.179769638 -0.05935534
## NUM_MOTORIST_KILLED                          -1.131633387 -0.88772171
## VEHICLE_14.DR.SEDAN                          -0.411465476 -0.30861113
## VEHICLE_1BOX.TRUCK                            0.926695635  0.10902693
## VEHICLE_1BUS                                  0.558998156  1.00527365
## VEHICLE_1OTROS                                0.641612799  0.60050456
## VEHICLE_1PASSENGER.VEHICLE                   -0.147110463 -0.17605129
## VEHICLE_1PICK.UP.TRUCK                        0.013877005 -0.27192311
## VEHICLE_1SEDAN                               -0.398715088  0.17213233
## VEHICLE_1SPORT.UTILITY...STATION.WAGON       -0.124820080 -0.57734601
## VEHICLE_1STATION.WAGON.SPORT.UTILITY.VEHICLE -0.410308042 -0.03140306
## VEHICLE_1TAXI                                 2.406097604 -0.60523382
## VEHICLE_1VAN                                  1.192997490  0.07636370
## VEHICLE_2                                    -0.310648489  0.31906128
## VEHICLE_2BIKE                                 0.751452567  0.16523582
## VEHICLE_2BUS                                  0.545434953 -0.41422432
## VEHICLE_2OTROS                                0.580756500  0.46476348
## VEHICLE_2PASSENGER.VEHICLE                   -0.231153104 -0.36552747
## VEHICLE_2PICK.UP.TRUCK                        0.164372904 -0.49973998
## VEHICLE_2SEDAN                               -0.254682722  0.38339403
## VEHICLE_2SPORT.UTILITY...STATION.WAGON       -0.062881235 -0.51537006
## VEHICLE_2STATION.WAGON.SPORT.UTILITY.VEHICLE -0.430921613 -0.30985802
## VEHICLE_2TAXI                                 2.732860997 -1.27570523
## VEHICLE_2UNKNOWN                             -0.135614769  0.37706043
## CAUSEBACKING.UNSAFELY                        -0.291195525  0.17418219
## CAUSEDRIVER.INATTENTION.DISTRACTION           0.083965448 -0.63961134
## CAUSEFAILURE.TO.YIELD.RIGHT.OF.WAY           -0.506731811 -0.99579476
## CAUSEFATIGUED.DROWSY                          0.288201981 -0.19498850
## CAUSEFOLLOWING.TOO.CLOSELY                   -0.005014399 -0.07208940
## CAUSEOTHER.VEHICULAR                          0.983710203  2.13525159
## CAUSEOTROS                                   -0.059111918  0.41208734
## CAUSEPASSING.OR.LANE.USAGE.IMPROPER           0.093391776 -0.74877957
## CAUSEPASSING.TOO.CLOSELY                      0.157325716  1.29819991
## CAUSETRAFFIC.CONTROL.DISREGARDED             -0.562986448  0.05403137
## CAUSETURNING.IMPROPERLY                       0.262584857  0.75013756
##                                                      LD3         LD4
## NUM_PERSONS_INJURED                          -0.00895853  0.02838755
## NUM_PERSONS_KILLED                           -1.50634719  0.92500320
## NUM_PEDESTRIANS_INJURED                      -0.76878255 -0.72717844
## NUM_PEDESTRIANS_KILLED                       -3.91548035 -3.17652726
## NUM_CYCLIST_INJURED                          -1.74058530 -0.59396703
## NUM_CYCLIST_KILLED                           -9.96263010  3.67370388
## NUM_MOTORIST_INJURED                          0.24394897  0.18831633
## NUM_MOTORIST_KILLED                           2.41828558  3.31300347
## VEHICLE_14.DR.SEDAN                           0.46152331 -0.60810832
## VEHICLE_1BOX.TRUCK                           -1.80365434 -0.04227749
## VEHICLE_1BUS                                 -0.45153957  1.01994454
## VEHICLE_1OTROS                               -0.52243962 -0.38113243
## VEHICLE_1PASSENGER.VEHICLE                   -0.16263732  0.02022199
## VEHICLE_1PICK.UP.TRUCK                        0.21177877  1.77948931
## VEHICLE_1SEDAN                                0.10510907  0.59289620
## VEHICLE_1SPORT.UTILITY...STATION.WAGON        0.34357441 -0.02493295
## VEHICLE_1STATION.WAGON.SPORT.UTILITY.VEHICLE  0.10647844 -0.82720416
## VEHICLE_1TAXI                                 0.84395809 -0.27800398
## VEHICLE_1VAN                                 -1.05193079  1.38741755
## VEHICLE_2                                     0.56272629  0.80601929
## VEHICLE_2BIKE                                -0.90235161  0.84625546
## VEHICLE_2BUS                                 -0.54307477 -0.12715813
## VEHICLE_2OTROS                               -0.87814306 -0.07157026
## VEHICLE_2PASSENGER.VEHICLE                    0.19633924 -0.42023699
## VEHICLE_2PICK.UP.TRUCK                        0.80101500  0.14084964
## VEHICLE_2SEDAN                                0.30512081 -0.21954382
## VEHICLE_2SPORT.UTILITY...STATION.WAGON        0.02175711  0.07516290
## VEHICLE_2STATION.WAGON.SPORT.UTILITY.VEHICLE -0.38749189 -0.14458771
## VEHICLE_2TAXI                                 0.11837183  0.08264591
## VEHICLE_2UNKNOWN                              0.02953143 -1.39827085
## CAUSEBACKING.UNSAFELY                        -0.11495147 -0.44339853
## CAUSEDRIVER.INATTENTION.DISTRACTION           0.20737994  0.25285191
## CAUSEFAILURE.TO.YIELD.RIGHT.OF.WAY           -0.12038936 -0.04490143
## CAUSEFATIGUED.DROWSY                         -1.37311710  1.56271121
## CAUSEFOLLOWING.TOO.CLOSELY                   -0.66752894  0.07949691
## CAUSEOTHER.VEHICULAR                          1.53472807 -0.30980607
## CAUSEOTROS                                    0.02437761  0.59390473
## CAUSEPASSING.OR.LANE.USAGE.IMPROPER          -0.65250182 -1.28969408
## CAUSEPASSING.TOO.CLOSELY                     -1.38202580 -0.27406315
## CAUSETRAFFIC.CONTROL.DISREGARDED             -0.13068444 -1.86290398
## CAUSETURNING.IMPROPERLY                       1.37009487 -1.18252349
## 
## Proportion of trace:
##    LD1    LD2    LD3    LD4 
## 0.8442 0.0797 0.0547 0.0213

4.3.11 5.2 Evaluamos el modelo:

predictions <- predict(lda_model, newdata = test_data)

# Matriz de confusión
table(Predicho = predictions$class, Real = test_data$BOROUGH)

##                Real
## Predicho        BRONX BROOKLYN MANHATTAN QUEENS STATEN ISLAND
##   BRONX            40       47        13     25             6
##   BROOKLYN        303      757       454    608            78
##   MANHATTAN       127      176       576    139            17
##   QUEENS          353      790       424    922           143
##   STATEN ISLAND     0        0         0      0             0

4.3.12 6.2 Visualizar las dos primeras funciones discriminantes:

lda_df <- data.frame(predictions$x, BOROUGH = test_data$BOROUGH)

ggplot(lda_df, aes(x = LD1, y = LD2, color = BOROUGH)) +
  geom_point(alpha = 0.5) +
  labs(title = "Proyección LDA: Accidentes por distrito sin LATITUDE y LONGITUDE",
       x = "Función discriminante 1",
       y = "Función discriminante 2") +
  theme_nyc()

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

Una vez eliminadas las variables LONGITUDE y LATITUDE del análisis, se aprecia una pérdida notable en la calidad de la agrupación, los clústeres resultantes no muestran una estructura tan clara como antes. Esto indica que los grupos formados tienen menos coherencia interna y mayor solapamiento entre ellos. En otras palabras, sin la información geográfica precisa que aportaban esas variables, el modelo tiene más dificultades para diferenciar patrones consistentes en los datos. Este resultado confirma que LONGITUDE y LATITUDE eran variables muy influyentes —quizá demasiado— y que, al eliminarlas, se muestra la verdadera dificultad del problema de segmentación.

4.4 Análisis de correspondencia

Para llevar a cabo el análisis de correspondencias, consideramos interesante estudiar si existía alguna relación entre la hora y el día en que ocurren los accidentes y el resto de variables. Para ello, fue necesario calcular e incorporar al conjunto de datos las variables HOUR y DAY_OF_WEEK:

data_sampled <- data_sampled |> 
  mutate(HOUR = as.numeric(format(strptime(TIME, format = "%H:%M"), "%H")),
         DAY_OF_WEEK = weekdays(as.Date(DATE, format = "%Y-%m-%d")))
data_sampled$HOUR <- as.factor(data_sampled$HOUR)
data_sampled$DAY_OF_WEEK <- toupper(data_sampled$DAY_OF_WEEK)
data_sampled$DAY_OF_WEEK <- factor(trimws(data_sampled$DAY_OF_WEEK))  # Eliminar los posibles espacios

4.4.1 1. Tipo de vehiculo y causa:

# Crear y limpiar la tabla de contingencia
contingency_table <- table(data_sampled$VEHICLE_1, data_sampled$CAUSE)
contingency_table <- contingency_table[rowSums(contingency_table) > 0, colSums(contingency_table) > 0]

# Realizar el análisis de correspondencia
ca_result <- FactoMineR::CA(contingency_table, graph = FALSE)

# Extraer coordenadas
row_coord <- as.data.frame(ca_result$row$coord)
col_coord <- as.data.frame(ca_result$col$coord)

# Añadir etiquetas y tipo
row_coord$label <- rownames(row_coord)
col_coord$label <- rownames(col_coord)
row_coord$type <- "Fila"
col_coord$type <- "Columna"

# Combinar datos
biplot_data <- rbind(
  data.frame(Dim1 = row_coord[,1], Dim2 = row_coord[,2], label = row_coord$label, type = row_coord$type),
  data.frame(Dim1 = col_coord[,1], Dim2 = col_coord[,2], label = col_coord$label, type = col_coord$type)
)

# Gráfico interactivo con Plotly
plot_ly(
  data = biplot_data,
  x = ~Dim1, y = ~Dim2,
  type = 'scatter',
  mode = 'markers+text',
  text = ~label,
  textposition = 'bottom center',
  marker = list(size = 7),
  color = ~type,
  colors = c("Fila" = "#1F3B73", "Columna" = "#FF4C4C"),
  textfont = list(size = 12)
) |> 
  layout(
    title = list(
      text = "<b>Biplot del Análisis de Correspondencias<b>",
      x = 0.5,
      font = list(size = 16, color = "#1F3B73", family = "Roboto Condensed")
    ),
    xaxis = list(
      title = list(text = "<b>Dim 1<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    yaxis = list(
      title = list(text = "<b>Dim 2<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    plot_bgcolor = "#F4F4F4",
    paper_bgcolor = "#F4F4F4",
    legend = list(
      bgcolor = "#F4F4F4",
      bordercolor = "#F4F4F4",
      font = list(color = "#2E2E2E", family = "Roboto Condensed")
    )
  )

4.4.2 2. Tipo de vehiculo y hora:

# Crear una tabla de contingencia
contingency_table <- table(data_sampled$VEHICLE_1, data_sampled$HOUR)

contingency_table <- contingency_table[rowSums(contingency_table) > 0, colSums(contingency_table) > 0]

# Realizar el análisis de correspondencia
ca_result <- FactoMineR::CA(contingency_table, graph = FALSE)

# Extraer coordenadas
row_coord <- as.data.frame(ca_result$row$coord)
col_coord <- as.data.frame(ca_result$col$coord)

# Añadir etiquetas y tipo
row_coord$label <- rownames(row_coord)
col_coord$label <- rownames(col_coord)
row_coord$type <- "Fila"
col_coord$type <- "Columna"

# Combinar datos
biplot_data <- rbind(
  data.frame(Dim1 = row_coord[,1], Dim2 = row_coord[,2], label = row_coord$label, type = row_coord$type),
  data.frame(Dim1 = col_coord[,1], Dim2 = col_coord[,2], label = col_coord$label, type = col_coord$type)
)

# Gráfico interactivo con Plotly
plot_ly(
  data = biplot_data,
  x = ~Dim1, y = ~Dim2,
  type = 'scatter',
  mode = 'markers+text',
  text = ~label,
  textposition = 'bottom center',
  marker = list(size = 7),
  color = ~type,
  colors = c("Fila" = "#1F3B73", "Columna" = "#FF4C4C"),
  textfont = list(size = 12)
) |> 
  layout(
    title = list(
      text = "<b>Biplot del Análisis de Correspondencias<b>",
      x = 0.5,
      font = list(size = 16, color = "#1F3B73", family = "Roboto Condensed")
    ),
    xaxis = list(
      title = list(text = "<b>Dim 1<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    yaxis = list(
      title = list(text = "<b>Dim 2<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    plot_bgcolor = "#F4F4F4",
    paper_bgcolor = "#F4F4F4",
    legend = list(
      bgcolor = "#F4F4F4",
      bordercolor = "#F4F4F4",
      font = list(color = "#2E2E2E", family = "Roboto Condensed")
    )
  )

4.4.3 3. Causa del accidente y hora:

# Crear una tabla de contingencia
contingency_table <- table(data_sampled$CAUSE, data_sampled$HOUR)
contingency_table <- contingency_table[rowSums(contingency_table) > 0, colSums(contingency_table) > 0]

# Realizar el análisis de correspondencia
ca_result <- FactoMineR::CA(contingency_table, graph = FALSE)

# Extraer coordenadas
row_coord <- as.data.frame(ca_result$row$coord)
col_coord <- as.data.frame(ca_result$col$coord)

# Añadir etiquetas y tipo
row_coord$label <- rownames(row_coord)
col_coord$label <- rownames(col_coord)
row_coord$type <- "Fila"
col_coord$type <- "Columna"

# Combinar datos
biplot_data <- rbind(
  data.frame(Dim1 = row_coord[,1], Dim2 = row_coord[,2], label = row_coord$label, type = row_coord$type),
  data.frame(Dim1 = col_coord[,1], Dim2 = col_coord[,2], label = col_coord$label, type = col_coord$type)
)

# Gráfico interactivo con Plotly
plot_ly(
  data = biplot_data,
  x = ~Dim1, y = ~Dim2,
  type = 'scatter',
  mode = 'markers+text',
  text = ~label,
  textposition = 'bottom center',
  marker = list(size = 7),
  color = ~type,
  colors = c("Fila" = "#1F3B73", "Columna" = "#FF4C4C"),
  textfont = list(size = 12)
) |> 
  layout(
    title = list(
      text = "<b>Biplot del Análisis de Correspondencias<b>",
      x = 0.5,
      font = list(size = 16, color = "#1F3B73", family = "Roboto Condensed")
    ),
    xaxis = list(
      title = list(text = "<b>Dim 1<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    yaxis = list(
      title = list(text = "<b>Dim 2<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    plot_bgcolor = "#F4F4F4",
    paper_bgcolor = "#F4F4F4",
    legend = list(
      bgcolor = "#F4F4F4",
      bordercolor = "#F4F4F4",
      font = list(color = "#2E2E2E", family = "Roboto Condensed")
    )
  )

4.4.4 4. Causa del accidente y día de la semana:

# Crear una tabla de contingencia
contingency_table <- table(data_sampled$CAUSE, data_sampled$DAY_OF_WEEK)
contingency_table <- contingency_table[rowSums(contingency_table) > 0, colSums(contingency_table) > 0]

# Realizar el análisis de correspondencia
ca_result <- FactoMineR::CA(contingency_table, graph = FALSE)

# Extraer coordenadas
row_coord <- as.data.frame(ca_result$row$coord)
col_coord <- as.data.frame(ca_result$col$coord)

# Añadir etiquetas y tipo
row_coord$label <- rownames(row_coord)
col_coord$label <- rownames(col_coord)
row_coord$type <- "Fila"
col_coord$type <- "Columna"

# Combinar datos
biplot_data <- rbind(
  data.frame(Dim1 = row_coord[,1], Dim2 = row_coord[,2], label = row_coord$label, type = row_coord$type),
  data.frame(Dim1 = col_coord[,1], Dim2 = col_coord[,2], label = col_coord$label, type = col_coord$type)
)

# Gráfico interactivo con Plotly
plot_ly(
  data = biplot_data,
  x = ~Dim1, y = ~Dim2,
  type = 'scatter',
  mode = 'markers+text',
  text = ~label,
  textposition = 'bottom center',
  marker = list(size = 7),
  color = ~type,
  colors = c("Fila" = "#1F3B73", "Columna" = "#FF4C4C"),
  textfont = list(size = 12)
) |> 
  layout(
    title = list(
      text = "<b>Biplot del Análisis de Correspondencias<b>",
      x = 0.5,
      font = list(size = 16, color = "#1F3B73", family = "Roboto Condensed")
    ),
    xaxis = list(
      title = list(text = "<b>Dim 1<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    yaxis = list(
      title = list(text = "<b>Dim 2<b>", 
                   font = list(size = 12, color = "#2E2E2E", family = "Roboto Condensed")),
      tickfont = list(color = "#2E2E2E", family = "Roboto Condensed"),
      gridcolor = "#DADADA"
    ),
    plot_bgcolor = "#F4F4F4",
    paper_bgcolor = "#F4F4F4",
    legend = list(
      bgcolor = "#F4F4F4",
      bordercolor = "#F4F4F4",
      font = list(color = "#2E2E2E", family = "Roboto Condensed")
    )
  )

5 CONCLUSIONES